## "Public Preferences for International Law Compliance: 
##  Respecting Legal Obligations or Conforming to Common Practices?"
##
##
## Reproduce: 
##      - Figure SI 6 
##
## Input:  
##      - data/data_cleaned.csv
##
## Output: 
##      - fSI6.pdf
##
## Saki Kuzushima
## 01/06/2023


data <- read.csv(file="data/data_cleaned.csv")
source("code/functions_fSI6.R")

data$treaty_binary <- NA
data$treaty_binary <- ifelse(data$treaty == 1 | data$treaty == 2,
                             1, data$treaty_binary)
data$treaty_binary <- ifelse(data$treaty == 3 | data$treaty == 4 | data$treaty == 5,
                             0, data$treaty_binary)


item <- c("surname3", "whale1","hate1", "death1")
res <- expand.grid(m_idx=c(1, 0), item=item, stringsAsFactors=F)
res$lower <- 0
res$point <- 0
res$upper <- 0
res$p <- 0
res$se <- 0

for (i in seq(1, nrow(res))){
  item <- res[i, "item"]
  m_idx <- res[i, "m_idx"]
  res[i, c("lower", "point", "upper", "p", "se")] <- 
    diffmean_block(data=data, item=item, 
                   m="treaty_binary", m_idx=m_idx)
}

# Multiple testing correction
res$p_bh <- p.adjust(res$p, method="BH")

# adaptive shrinkage
ash_fit <- ash(betahat=res$point, 
               sebetahat=res$se, 
               prior="uniform", 
               mixcompdist="unif",
               grange=c(-5, 5),
               outputlevel=2)
ash_result <- ash_fit$result
ash_ci <- ashci(ash_fit)
res_ash <- data.frame(m_idx=res$m_idx,
                      item=res$item,
                      lower=ash_ci[,1], 
                      point=ash_result[,"PosteriorMean"],
                      upper=ash_ci[,2])
x_index_temp <- c(0.625, 0.875, 1.5, 1.5+0.625, 1.5+0.875, 3) 
x_index <- c(0, x_index_temp, 3+x_index_temp, 6+x_index_temp, 9+x_index_temp)

# fiugre 
point <- c(NA,res$point[1], res_ash$point[1], NA, 
                res$point[2], res_ash$point[2], NA, 
                res$point[3], res_ash$point[3], NA, 
                res$point[4], res_ash$point[4], NA,
                res$point[5], res_ash$point[5], NA, 
                res$point[6], res_ash$point[6], NA, 
                res$point[7], res_ash$point[7], NA,
                res$point[8], res_ash$point[8], NA)
upper <- c(NA, res$upper[1], res_ash$upper[1], NA, 
                res$upper[2], res_ash$upper[2], NA, 
                res$upper[3], res_ash$upper[3], NA, 
                res$upper[4], res_ash$upper[4], NA,
                res$upper[5], res_ash$upper[5], NA, 
                res$upper[6], res_ash$upper[6], NA, 
                res$upper[7], res_ash$upper[7], NA,
                res$upper[8], res_ash$upper[8], NA)
lower <- c(NA, res$lower[1], res_ash$lower[1], NA, 
                res$lower[2], res_ash$lower[2], NA, 
                res$lower[3], res_ash$lower[3], NA, 
                res$lower[4], res_ash$lower[4], NA,
                res$lower[5], res_ash$lower[5], NA, 
                res$lower[6], res_ash$lower[6], NA, 
                res$lower[7], res_ash$lower[7], NA,
                res$lower[8], res_ash$lower[8], NA)
point_style <- c(NA,1,16,NA,1,16,NA,2,17,NA,2,17,NA,0,15,NA,0,15,NA,5,18,NA,5,18,NA)

pdf("results/main_si/fSI6.pdf")
par(mfrow=c(3,1), mar=c(1,4,1,1)+0.1, oma=c(2,1,3,1)+0.1)
plot(x=x_index, y=point, type="n",
     xlim=c(0.5, 11.5), ylim=c(-0.5, 0.5),
     pch=point_style, cex=1.5,
     yaxt="n", xaxt="n", xlab="", ylab="Constitution - Int'l law ",
     cex.lab=1.25,
     col="gray50")
abline(h=0, col="red")
abline(v=3, lty="dotted")
abline(v=6, lty="dotted")
abline(v=9, lty="dotted")
axis(side=3, at=c(1.5, 1.5+3, 1.5+6, 1.5+9),
     labels=c("Surname", "Whale", "Hate Speech", "Death Penalty"),
     tick=F,  line=0, outer=T, cex.axis=1.25)
axis(side=3, at=c(0.75, 2.25, 0.75+3, 2.25+3, 0.75+6, 2.25+6, 0.75+9, 2.25+9),
     labels=rep(c("HI", "LO"),4),
     tick=F,  line=0, cex.axis=1.15)
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       y=point[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       pch=point_style[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       cex=c(rep(1.5,12), 1.5,2,1.5,2),
       col="grey50") 
segments(x0=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)],
         x1=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y0=lower[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y1=upper[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         cex=1.5,
         col="grey50") 
dev.off()



